A variational pseudo-self-interaction correction approach: ab-initio description of 

correlated oxides and molecules 
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We present a fully variational generalization of the pseudo self-interaction correction (VPSIC) 
approach previously presented in two implementations based on plane-waves and atomic orbital 
basis set, known as PSIC and ASIC, respectively. The new method is essentially equivalent to the 
previous version for what concern the electronic properties, but it can be exploited to calculate 
total-energy derived properties as well, such as forces and structural optimization. We apply the 
method to a variety of test cases including both non-magnetic and magnetic correlated oxides and 
molecules, showing a generally good accuracy in the description of both structural and electronic 
properties. 
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I. INTRODUCTION 

The ab-initio, density-functional theory (DFT) based 
determination of structural and electronic properties of 
correlated systems remains up today an outstanding chal- 
lenge. The most relevant bottleneck rests in the fact 
that advanced approaches that are general and power- 
ful enough to tackle the description of strong-correlated 
systems are also usually heavily demanding in terms of 
required computing resources. This limits the system size 
that can be practically afforded to few tents of atoms per 
unit cell. On the other hand, most of the interesting 
behaviours encountered in correlated systems requires 
large cell size already at the level of bulk properties, 
due the possible coexistence of several juxtapposed order- 
ings (structural, magnetic, orbital and charge ordering). 
Furthermore, several much celebrated phenomenologies 
involve doping (high-Tc superconductivity in cuprates, 
colossal magnetoresistivity in manganites, magnetic or- 
dering in diluted semiconductors, etc.) whose treatment 
at generic concentration is a formidable computing task; 
finally, oxide interfaces and multilayers which are ground 
of recent intriguing discoveries such as 2D electron gas 
behaviour 1 may require even larger simulation effort. 

As a matter of fact, we are in the need of treating 
about one or two hundred atom systems, a size that can 
be only afforded at a computational cost substantially 
similar to that of standard local (spin) density functional 
(L(S)DA) theory or its generalized gradient (GGA) ver- 
sion. Bcyond-LSDA approaches which are agile enough 
to satisfy this requirement are very few: the very popu- 
lar LDA+U^ is certainly one of those; a later approach 
with similar characteristics is the PSIC, implemented in 
the past in two different setting: plane-wave and ultra- 
soft pseudopotential basis se t 49 ' 50 in the framework of the 
home-made PWSIC code; local orbital and pseudopo- 
tential basis set (also called ASIC)^ 1 - in the framework 
of the SIESTA code. LDA+U and PSIC/ ASIC move 
from different conceptual viewpoints: the former intro- 



duces an effective Coulomb repulsion which is tipically 
mistreated in LSDA; the latter subtracts from the LSDA 
functional an approximate (i.e. atomic orbital-averaged) 
self-interaction (SI), that is the unphysical interaction of 
an electron with its own generated potential. 

Despite this difference, the two theories act in similar 
way, i.e. correcting the LSDA eigenvalues by a quan- 
tity linearly dependent on the orbital occupations; in 
fact, the PSIC can be substantially viewed as an all- 
orbital generalization of the LDA+U, but with two rele- 
vant advantages over the latter: the capability of cur- 
ing the LSDA failures in more general situations (i.e. 
not limited to magnetic and insulating systems) and 
the absence of an explicit parametric dependence. In 
PSIC/ASIC, indeed, the function of U si replaced by the 
atomic self-interaction potentials, which are extracted 
from the free atom (thus they are universal, i.e. only 
dependent on the atomic species of the systems) and 
then plug into the band structure Hamiltonian. At vari- 
ance with more fundamental SI removal strategies such 
as the original Perdew-Zunger approach 64 (PZ-SIC), or 
the later generalization to extended system o 59 ' 60 which 
draw on dramatic complications in formalism and con- 
ceptual interpretations^^, the PSIC keep all the sim- 
plicity tipical of LDA/GGA: a single-particle potential 
which is orbital-independent and translationally invari- 
ant (i.e. obeying Bloch symmetry), and energy func- 
tional invariant under unitary rotation of the occupied 
eigenstates. 

The PSIC/ASIC approach demonstrated a consistent 
accuracy in the description of the electronic properties 
for a vast range of systems^ -, at a computational cost 
not much larger than that of the LSDA itself (in particu- 
lar the ASIC is the implementation of choice for large-size 
systems as it can easily afford few hundreds atom simula- 
tions, while the heavier PSIC can be seen as the standard 
reference for what concern the methodological accuracy). 
Despite these virtues, PSIC and ASIC have had a rela- 
tively small following with respect to the LDA+U so far, 
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mainly in reason of the serious hindrance represented by 
the lack of variationality: in Refl49lthe PSIC potential 
is in fact generated as an ansatz on the single-particle 
Kohn-Sham (KS) potential, not from a germinal energy 
functional. This shortcoming precludes, in practice, the 
access to ground-state total energy properties. 

In this paper we point to overcome this limitation, pre- 
senting a fully variational version of PSIC/ASIC (here- 
after called VPSIC, to be distinguished by their pre- 
vious siblings) built out of an energy functional which 
recovers, by Euler-Lagrange derivative, a set of single- 
particle KS equations essentially similar to those of the 
elder PSIC/ASIC scheme, thus keeping the successful de- 
scription of the electronic properties, but adding up the 
possibility to deliver all those ground-state properties ex- 
pected by standard ab-initio theories. 

Here we describe the new formulations and give wide 
evidence of the aforementioned capabilities presenting re- 
sults for a range of bulks (including non magnetic and 
magnetic titanates and manganites) and molecules. Ex- 
tended systems are treated using the plane-wave basis 
implementation, the molecules by the atomic orbital ba- 
sis implementation. We remark that other VPSIC results 
have been already presented in separate recent publica- 
tions: a lenghtly description of the properties of transi- 
tion metal monoxides^; a detailed account of the proper- 
ties of 2D electron gas formation at the SrTiOs/LaAlOa 
interface^; together with the present article, they furnish 
a solid evidence of the effectiveness of the new method in 
the study of generic systems characterized by strong to 
moderate electron correlation. 

The article is organized as it follows: in Section HIl the 
general formulation is illustrated; in Section Mil results 
for non magnetic oxide insulators (Subsection 1111 Al) . ti- 
tanates (Subsection IIII Bl) and manganites (Subsection 
IHICj) are presented. Section |TV] is devoted to illustrate 
results for molecules. Finally, in Section [V] we draw our 
summary and conclusions. Several specific implementa- 
tions are included in the Appendices: in I VII and IVIII we 
present an extension of the VPSIC energy functional and 
forces formulations, respectively, for the case of ultrasoft 
pseudopotentials (USPP), while IVnTl and ?? are dedi- 
cated to the description of VPSIC energy functional and 
forces formulation specific for atomic orbital basis set. 



II. VARIATIONAL 
PSEUDO-SELF-INTERACTION FORMULATION 

In this Section we present the generic variational for- 
mulation, not related to any specific basis function im- 
plementation. 



A. VPSIC energy functional and related 
Kohn-Sham equations 

We start from the following VPSIC energy functional: 



E vpsic m] = E Ls °m] t e s§i„ P ° iu (i) 

where E LSD is the usual LSDA energy functional: 

E LSD [{^}]=T s [{^}]+E H [n]+E xc [n + ,n_]+E ton [{4,}] 

written as sum of (non-interacting) kinetic (T s ), 
Hartree (E#), exchange-correlation (E xc ), and electron- 
ion (Ei ora ) energies (ip are single-particle wavefunctions, 
n+ and n_ up- and down-polarized electron densities, 
and n=n + +n_). EqQ] follows the spirit of the original 
Perdew-Zunger procedure^ (hereafter called PZ-SIC), 
and subtracts from the LSDA total energy a SI term 
written as a sum of effective single-particle SI energies 
(£ SI ) rescaled by some orbital occupations P. Here i, j 
are sets (li,nii, U^rrij) of atomic quantum numbers (typ- 
ically relative to a minimal atomic wavefunction basis 
set) while a and v indicate spin and atomic site, respec- 
tively (non-diagonal formulation is necessary to enforce 
covariancy, thus i={li 1 rrii), j={li,m,j)). 

Most of the peculiarity of the VPSIC approach resides 
in the way the second term of Eqfl] is written for an 
extended system whose eigenfunctions are Bloch states 
(Vnk)- The orbital occupations are then calculated as 
projection of Bloch states onto localized (atomic) orbitals 
(hereafter indicated as {</>}): 

Piju = E /nk «kl^V) to>|<k>, (2) 
nk 

where fZ± are Fermi occupancies. For the effective SI 
energies we adopt a similar approach: 

£ i/«u = E fSk WSkhi,»)Gj ( T^Kk) (3) 

nk 

where 7^,, is the projection function associated to the 
SI potential of the i th atomic orbital centered on atom v. 

•y iv {r-R v ) = VH XC [niu{r-~Rv),Q}<t>iv{r-'R v ), (4) 

where n,i V (r) =4>% (r). We express the Hartree plus 
exchange-correlation atomic SI potential Vhxc m ra- 
dial approximation: V Hxc = V H [ni iU ] + V xc [ni iU , 0] = 
dEHxc[ni iV \/dni iV (calculated at full polarization: n=n + , 
n_=0). Finally, the Cij are normalization coefficients: 

C^ mumj = |dr^ mi (r)VHxcK,,(r),0]^ mj (r) (5) 

with h=lj. They are purely atomic and do not de- 
pend on the atomic positions. The use of projectors 7 in 
Eq|3]is aimed at casting the SI energy in fully non-local 
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form (analogous to the fully-non local pseudopotential 
form due to Kleinman and Bylander— ) which consent a 
huge saving of computational effort when calculated in 
reciprocal space. 

To grab the idea behind Equations |31 IH and [SJ no- 
tice that in the limit of large atomic separation, the 
Bloch states V'nk are brought back to atomic orbitals <f>i V , 
and £>fj av to atomic SI energies ef 1 (discarding spin and 
atomic index): 



rSl 



drmiv) (VffK(r)] +14cK(r),0] 



(6) 



Thus, the orbital occupations Pg- V (if suitably normal- 
ized) act as scaling factors for the atomic SI energies, 
assumed as the upper limit of the SI correction ampli- 
tude. 

We remark that in the atomic limit EqUgoes back to 
the PZ-SIC total energy espression only for what concern 
the Hartree SI part, while our SI exchange-correlation 
energy density ((l/2)V xc [ni, 0]) departs from the PZ-SIC 
espression e xc [rii,0], since V xc = e xc + mde xc /drii). 

From Eq[T] we obtain the corresponding VPSIC KS 
Equations through the usual Eulcr-Lagrange derivative: 



spectrum substantially similar to that of the non varia- 
tional scheme, but with the added bonus to derive from 
the VPSIC energy functional. 

In DFT methods it is customary to rewrite the total 
energy in terms of eigenvalue sum. Indicating with e n ka- 
the LSDA eigenvalues, it is immediate to verify that: 



/nk ^nko 

nkcr 



E /"k «l 



nker 



dE VPSIC s 



then Eq. [T]can be rewritten as: 



(11) 



E v™c m] = ^[{ffl + \ J2 e&y Pliu (12) 



where: 



E LbU l{iP}] = > f^inka +E Hxc [n+(r),n-(r)} 



dE 



VPSIC 



ka 



dE, 



si 



dPZ, 



dtb* 

rnka 



where e n fc CT are VPSIC eigenvalues, and; 



l Lsn { 



r = 



+ V H [n(r)} + V xc [n+{r),n-(r)} + V ion (r) 
is the usual KS LSDA Hamiltonian, and: 



d£ SI 



(8) 



(9) 



E„ 



drn CT (r)^ c K(r),n_(r)] (13) 



is the LSDA energy functional but now including the 
VPSIC eigenvalues in place of the LSDA eigenvalues. Fi- 
nally, in the original PSIC formulation the SI Vh X c poten- 
tial is rescaled by a relaxation factor a— 1/2, to take into 
account the screening (i.e. the suppression) of atomic 
self-interaction caused by the surrounding charge of the 
other atoms (see Ref.— for an extended discussion). A 
careful testing on a large series of compounds^ confirmed 
that this relaxation value is the most adequate for almost 
all the examined bulk systems, whereas for molecules the 
atomic SI (a=l) is the most appropriate choice. We man- 
tain this recipe in the present formulation, using the 1/2 
factor only for calculations related to extended systems. 



B. Simplified variants of VPSIC and relation with 
the original non-variational method 



— ^ = |&,„) (^Kk). (io) 

Vnka 

The first sum term in curl parenthesis in EqJT] corre- 
sponds to the SI potential projector written as in the 
original VPSIC KS Equations. Since the two sums in 
curl parhenthesis are two ways to describe substantially 
the same physical quantity (i.e. the SI potential), in 
practice they give similar results when applied onto the 
Bloch state. It follows that Eq|7| describes an energy 



From the general espression of EqQ] two interesting 
subcases can be obtained: assumig fixed (i.e. non self- 
consistent) orbital occupations Py, in EqJT] the second 
term in curl brackets vanishes and the VPSIC-KS equa- 
tions reduce to those of the original PSIC scheme of 
Refl49l (indeed, it was previously pointed oub££ that the 
original scheme becomes variational at fixed orbital oc- 
cupations) . 

Another useful subcase is obtained replacing Eq|3]with 
a simplified expression: 
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3<T SI 



col _ pa i 



pa SI 



(14) 



where the atomic (in radial approximation) is given 
by EqEl Using EqQ~4| previous EqsfflandE] 

E VPSic 0[m = E LSD [m-\Y,P^ P ^^ as) 



projector centered on v. This is not necessarely true if the 
orbital occupations are to be re-orthonormalized on the 
cell. This choice, which complicates sensitively the above 
formulation, will be discussed in detail in the Appendices, 
together with the generalization of the method to either 
plane-waves plus USPP or local-orbital plus pseudopo- 
tential approach. 



III. RESULTS: EXTENDED SYSTEMS 



E 



LSD 



It. 



per pa SI 

ijis jiv jv 



Q E VPSIC 

dtb*. 



= h 



LSD 



^Pnka 



E 

ijv 



dP 



Jiv SI 
* jv 



(16) 



(17) 



These simplified VP SIC formalism (hereafter indicated 
as VPSICo) is a computationally convenient alternative 
(especially in terms of required memory) to perform 
structural optimizations in large size-systems. In the re- 
sults Section we will give evidence that at least in case of 
non-magnetic semiconductors and insulators it furnishes 
electronic and structural properties in good accord with 
the VPSIC. However it is tipically less satisfying for the 
electronic properties of magnetic systems. 



C. Forces formulation 

In VPSIC the atomic forces formulation follows from 
the usual Hellmann-Feynmann procedure. It is obtained 
as the LSDA expression augmented by a further additive 
contribution due to the atomic-site dependence of the SI 
projectors: 



dE ypsic [W] 



? LSD 



~ E /:k{<Cki^)^-<7,>i< k >^[W]+c.c.} 

ij,nk(T ' 



^ frik 
ij,nk<r 



8K U 



X<MCk)+cc fl8) 



whereas in the simplified version, we have, in addition 
to F^ SD the quantity: 



E fn 



ij.nkcr 



SI 



BR. 



: )(^>Kk) + " (19) 



In writing Eqsfl8] and HH we have assumed that the 
force on a given atom v only depends on the single atomic 



We have pointed out in Section |TT] that there is a sub- 
stantial formal similarity between the KS Equations de- 
rived for the VPSIC and the elder PSIC/ASIC KS equa- 
tions. Our tests assess that the long series of results for 
the electronic properties obtained with the latter in the 
last few years remain valid even in the framework of the 
new theory, which gives small differences in e.g. band en- 
ergies and density of states. Hereafter we consider as test 
cases for VPSIC, materials either never tackled before 
(it is the case of titanates), or afforded in the past (e.g. 
CaMnOa) but now revisited to specifically address total 
energy-derived properties, such as equilibrium structure 
and magnetic exchange-interactios. 

Specifically, we selected three classes of systems which 
well represent the broad spectrum of the VPSIC capabil- 
ity, and at the same time are interesting compounds in 
itself either from conceptual or technological viewpoint: 
wide-gap oxide insulators, magnetic titanates representa- 
tive of 3d t2 g Mott-insulating perovskites, and magnetic 
manganites representative of 3d e g charge-transfer insu- 
lating perovskites. 

Regarding our technicalities, calculations are carried 
out with ultrasoft pseudopotentials^ and a plane-wave 
basis set with cut off ranging from 30 to 35 Ryd de- 
pending on the specific system, 6x6x6 special k-point 
grids for self-consistent calculations, 10x10x10 special 
k-points and linear tetrahedron interpolation method for 
density of states. The Ceperley-Alder-Perdew-Zunger 
local-density approximation is used for the exchange- 
correlation functional. Structural relaxations are carried 
out with a convergency threshold of 1 mRy/Bohr on the 
calculated forces. 



A. Wide-gap insulators: LaA10 3 , SrTi0 3 , Ti0 2 

As prototypes of non-magnetic wide-gap insulators we 
selected TiC>2 rutile and two perovskite oxides, LaA103 
and SrTiOs in their cubic bulk structure. Nowadays 
these materials are ravely investigate for their potential 
use in the field of functional oxides, thus the abundance 
of theoretical and experimental results make them ideal 
test cases for innovative theoretical methods. 

For non magnetic oxides the level of accuracy of stan- 
dard LDA calculations may vary according to the ex- 
amined property: typically a good rendition of struc- 
tural properties is juxtaposed to an unsatisfactory match 
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TABLE I: Direct (AE d ) and indirect (AEi) band gap en- 
ergies and O 2p manifold valence bandwidth (Wop) in eV, 
calculated within LDA, PSIC, and PSICo, compared to the 
experimental values. 



LDA 



PSIC 



PSIC, 







LDA PSIC PSICo 


Expt. 


Ti0 2 


AEi (r-L) 


1.88 


2.90 


2.59 






AE d (r) 


1.84 


2.93 


2.62 


[3.0,3.1^ 




Wop 


6.0 


6.5 


5.7 


~ t£ 


SrTiOs 


AEi (M-r) 


1.69 


2.94 


2.62 


3.2C 1 ' 




AEd (r) 


2.04 


3.30 


2.95 


3.75^ 




Wop 


5.0 


5.5 


4.8 


~ 6^ 


LaA10 3 


AEi (M-r) 


2.83 


5.23 


4.61 






AE d (r) 


3.17 


5.51 


4.89 


5.6 




Wop 


7.6 


8.38 


7.27 


~ 8-9^ 



of the calculated band structure and interband transi- 
tion energies, involving the well-known underestimation 
of the fundamental band gap, and the poor description 
of transition-metal d and (on a minor extent) oxygen 
p states. Our results will demonstrate that the VP- 
SIC achieves a sistematical improvement of the electronic 
properties over LSDA, while preserving a substantially 
similar accuracy for what concerns structural properties. 
However, it must be kept in mind that the detail of the re- 
sults may (and usually does) depend crucially on a num- 
ber of computational technicalities typical of our super- 
cell simulations, which goes beyond the formalism (e.g. 
the type of wavef unction basis set, the type of used pseu- 
dopotentials, etc.). Thus, in order to reach an unbiased 
evaluation it is always useful to refer the PSIC results 
not only to the experiment values but even to their LDA 
counterpart, obtained in condition of identical technical- 
ities. 

We start with the analysis of the band energies of TiC>2 
in FiglU where results for VPSIC, VPSIC , and LDA are 
reported (for an unbiased comparison of the methods we 
consider the bands calculated at the same (experimental) 
structure). The corresponding band gap energy values 
are reported in Tab|H As expected, TiC>2 shows the en- 
ergy gap between occupied O 2p valence and empty Ti 
3d conduction bands. According to VPSIC, the minimum 
gap is indirect between T at valence band top (VBT) and 
M (i.e. the BZ edge along [110]) at the conduction band 
bottom (CBB), while the direct gap is at T. The energy 
gap values are in satisfying agreement with the exper- 
iments, while the LDA results present the well-known 
band gap underestimation of nearly 40%, a typical LDA 
error bar for non-magnetic insulators. The PSICo, on the 
other hand, operates a partial (about 80%) recovery of 
the correct energy gap over the LDA result. While this 
may be somewhat unsatisfying for the purpose of predict- 
ing accurate photoemission energies, it may be sufficient 
for performing structural optimization at a pace substan- 
tially similar to that of the LDA itself. Notice also that 
for what concern the (mainly O 2p) valence bandwidth, 
VPSIC and VPSICo changes the LDA value (~ 6.0 eV) 
in opposite directions (~ 6.5 eV for VPSIC, 5.7 eV for 
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FIG. 1: Band energies of Ti02 (rutile) at experimental struc- 
ture (a=4.59 A, c/a=0.664, u=0.305) calculated by LDA, VP- 
SIC and VPSICo. K-points coordinates (units of n/a—n/b, 
tt/c) are r=(0,0,0), M=(1,1,0), R=(l,l,l). 



VPSICo); this illustrates the fact that the VPSIC is not 
simply a rescaled VPSIC. 

It is worthy emphasizing that the capability of VP- 
SIC (or VPSICo) to correct the LDA failure in case of 
non-magnetic oxides whose ground state is only residu- 
ally affected by d-type orbitals, spurs from its all-orbital 
corrective character (in particular from the O 2p band 
correction). In contrast, the GGA+U band gap of 2.2 
eV (with U=3.4 eV applied to Ti 3d) calculated in RefH 
is only marginally larger than the respective GGA value 
of 1.9 eV. The necessity of applying a corrective U onto 
the O 2p orbital energies was indeed discussed in previous 
works^. 

The band energies for bulk cubic SrTiOa is reported in 
Figs f51 The dominant (and eventually the secondary) or- 
bital character for each group of bands is also reported in 
Figure. The energy gap opens between CBB of mainly 
Ti 3d character, with higher contributions from Sr 4d 
states, and VBT bands dominated by O 2p states. Wee 
can see in Tab [I] that the LDA band gap is only ~ 55% of 
the experimental direct (3.25 eV) and indirect (3.75 eV) 
gap. The VPSIC, on the other hand, recovers most (~ 
90%) but not all the LDA discrepancy. This is in part 
attributable, rather than to a VPSIC insufficiency, to a 
much too low LDA value, due to our specific technical- 
ities and used pseudopotentials. Indeed previous LDA 
determinations^ gave 1.9 eV and 2.24 eV for indirect 
and direct band gap.and a "scissor" operator of 1.5 eV 
was employed in Ref|5j to readjust band energies with el- 
lipsometry data. According to our calculations the VP- 
SIC operate a "scissor" of about 1.3 eV with respect to 
the LDA band gap. Notice also that while in general 
the action of VPSIC over the LDA bands may be larg- 
erly k-point dependent, for these wide-gap, highly ionic, 
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r m r rr m r rr m r r 

FIG. 2: Band energies of cubic SrTi03 at experimental struc- 
ture (a=3.92 A). K-points coordinates (units of ir/a) are 

r=(o,o,o), m=(i,i,o), R=(i,i,i)- 



non-magnetic oxides the LDA band dispersion is not dra- 
matically modified, indeed, and the concept of scissor op- 
erator can be grossly justified "a posteriori". However, 
if the shape is not much changed, significant differences 
are visible in the bandwith: consider e.g. the first unoc- 
cupied doublet of Ti 3d t2 g character: in LDA it spans 
about 2.5 eV, while in VPSIC or VPSIC it is reduced to 
about 2 eV, in consequence of the enhanced charge local- 
ization (i.e. reduced p-d hybridization) which tipically 
follows from the removal of SI. On the other hand, the 
VPSIC increases the LDA O 2p bandwidth as a conse- 
quence of the different spectral weight distribution (the 
top manifold is purely O 2p, while the bottom (R point) 
shows a significant mixture of Ti 3d and Sr 3s states. 
The different effective SI energies related to these differ- 
ent orbital characters eventually stretches the band bot- 
tom down to lower energies, thus resulting in an increase 
of bandwidth. In VPSICo on the other hand the effective 
SI energies are fixed to their atomic counterparts, thus 
the effect over the LDA values is more similar to that of 
a rigid band shift. 

Finally, in Fig. [3]the band energies of cubic LaA103 are 
reported. Now the energy gap is between O 2p valence 
bands and mainly Al s,p conduction bands. The absence 
of 3d states produces an observed gap (5.6 eV) much 
larger than in SrTiOs, and again grossly underestimated 
in LDA (^56% of the experimental value). This case is 
prototypical in demonstrating the necessity to repair the 
LDA unrespectively on the presence of 3d states. Again, 
VPSIC and VPSICo works quite nicely in recovering a 
satisfying energy gap. However, as in case of SrTi03 they 
act in different ways for what concern the band width: 
VPSIC stretches the bottom of valence O 2p manifold 
down to lower energies, while in VPSICo the same bands 
span an energy window even smaller than in LDA. For 
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FIG. 3: Band energies of cubic LaAlOs at experimental struc- 
ture (a=3.82 A). K-points are the same as in Fig[2] 



TABLE II: Lattice parameters (in A) calculated within LDA, 
VPSIC, and VPSICo, compared to the experimental values. 
For SrTiOs and LaAlOs the cubic symmetry is assumed, while 
Ti02 is in tetragonal (rutile) structure. For T1O2 c/a is fixed 
to the experimental value 0.644, while the internal parameter 
u is calculated. 





LDA 


PSIC 


PSICo 


Expt. 


SrTiOa 


3.99 


3.97 


4.02 


3.92 


LaAlOs 


3.76 


3.74 


3.83 


3.82 


Ti0 2 a 


4.67 


4.62 


4.69 


4.59 


Ti0 2 u 


0.3021 


0.3066 


0.3066 


0.3048 



what concern the lower lying O s bands (between -16 eV 
and -18 eV) and La p (semicore) bands (around -14 eV), 
LDA and VPSIC give a quite similar description, while 
the two groups are overlapping according to VSICo- 

In Fig|?]we report total energies vs. lattice parameter 
for cubic SrTiOa, LaA103 and tetragonal Ti02, calcu- 
lated within LDA, VPSIC, and VPSIC for the three 
examined oxides. The calculated equilibrium parameters 
are reported in Tab|Hl in comparison with the experi- 
mental values. We can start the analysis considering the 
reference LDA values: as mentioned, LDA is tipically 
satisfying for what concern structural properties, but the 
level of accuracy may vary depending on technicalities. 
Thus while LDA is known to generally underestimate by 
1-2 % the equilibrium lattice constant, we have values 
for SrTiOs and T1O2 which overestimate the experiments 
by ~ 1.5%. On the other hand, the lattice constant of 
LaA103 is underestimated by little more than lsystem- 
atic (0.5-1%) reduction of the corresponding LDA values, 
thus ending up showing a very satifying accord with the 
experiment (within lgoes again in the opposite direction, 
as it increases sistematically the LDA values by 0.5%- 1%. 

The different behaviour of VPSIC and VPSIC can 
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LaA10 3 SrTi0 3 TiQ 2 




lattice constant (A) 

FIG. 4: (Color on-line) Total energies as a function of lattice 
parameter for for cubic SrTiOa, LaA103 and tetragonal Ti02, 
calculated by LDA, VPSIC, and VPSICo. 

be easily linked back to what we have commented for 
the bands: the VPSIC tend to stretches the occupied va- 
lence band manifold (O 2p, O 2s) with respect of the LDA 
values, thus reducing the effective screening and in turn 
the bond lenght. On the contrary, the VPSICo shrinks 
the occupied band manifolds with respect to the LDA, 
thus causing an increase of effective screening and bond 
length. The shrinking of equilibrium volume caused by 
the VPSIC was also reoprted for transition metal monox- 
ides MnO and NiO££. Clearly, this should not be in- 
tended as a 'universal' trend, as bandwidth and charge 
localization are case-dependent quantities, and as such 
are the VPSIC modifications over the LDA electronic 
properties. 

We can conclude this section emphasizing the overall 
good quality of the VPSIC predictions for the three ex- 
amined compounds, with lattice parameters and energy 
band gaps within 1-2% and ~10% from the experiments, 
respectively. While many more results will be necessary 
for a full assesment of the theory, what we have showed 
in this Section is sufficient to encourage the use of VPSIC 
for the investigation of wide-gap oxide insulators. 



B. Magnetic titanates: YTi0 3 , LaTi0 3 

Titanates characterized by the nominal Ti d 1 config- 
uration rank among the most peculiar and intriguing 
magnetic perovskites. At variance with the more inves- 
tigated classes of manganites and cuprates whose fun- 
damental chemistry is governed by 3d e 9 states, in ti- 
tanates the 3d valence states are 3d t2 9 , thus orbitals 



not directly oriented towards the oxygens; this produces 
a much weaker p-d hybridization and crystal field split- 
ting than in e g systems. However, experiments show that 
the phenomenology of these systems may be crucially af- 
fected by small structural details. 

A nice illustration of this over-sensitive magnetostruc- 
tural coupling is furnished by the compared study of 
YTi0 3 (YTO) and LaTi0 3 (LTO): both systems are 
Pnma perovskites, with relatively small Jahn- Teller (JT) 
distortions and large GdFeO-type octahedral rotations; 
the difference in cation size (with La bigger than Y) 
makes the amplitude of distortions and rotations sensibly 
wider in YTO than in LTO (in agreement with the well 
known space- filling criterion), in turn leading to quite 
a different magnetic behaviour: YTO is ferromagnetic 
(FM)i2riS with low Curie temperature T c =30 K, size- 
able band gap (~ 1 eV) and magnetic moments M=0.8 
/is, in line with a Ti d 1 ionic configuration; LTO, on 
the other hand, is antiferromag netic (AF) G-type 2 ^ with 
Tat=130 K, very small energy gap (~ 0.3 eV) and sen- 
sibly smaller magnetic moments (~ 0.57 //b)2£. A long- 
standing debate was ignited on the attempt to give an 
explaination to this much reduced magnetic moment and 
nearly isotropic spin- wave dispersion in LaTiOgi^. It was 
pointed out that a single electron in the triple-degenerate 
t% g manifold can fluctuate giving rise to an exotic 'orbital 
liquid' stat o 18 i 25 . However this fashinating hypothesis is 
contrasted by a series of evidence a 18 i 20 ~ 23 i 26 ' 29 that crys- 
tal field splitting is actually not small enough to keep 
the t2g degeneracy substantially unlifted, and instead a 
Jahn- Teller distorted, orbital-ordered state is realized in 
LTO, as well as in YTO. 

Needless to say, these issues stimulated a long series of 
attempts to describe the titanates by a variety of ab-initio 
approaches, including LSDA 30 , LDA+L T 27 ! 31 , and several 
LSDA+DMFT implementations^ 3 ^ 3 .. While our de- 
scription reproduces, at least in part, some of the previ- 
ous findings, one aspect makes our results very valuable: 
they capability to furnish a coherent rendition of struc- 
tural and electronic properties of on a purely ab-initio 
ground, in the framework of the same theory, and with- 
out inclusion of system-dependent parameters (e.g. U, 
J). 

We will proceed as following: we start illustrating the 
electronic properties of YTO, calculated at the experi- 
mental structure, as prototype of 'basic' t2 9 system, then 
we move to discuss the more peculiar LTO, highlighting 
the differences with respect to YTO, and finally we end 
up with the structural properties of both systems, which 
will furnish a rationale to their different behaviour. No- 
tice that we do not include any comparison with local- 
density results in the analysis: LSDA does not reproduce 
the Mott-insulating behaviour for these systems and in 
fact predict an unphysical non magnetic, metallic elec- 
tronic ground-state which is not meaningful for the pur- 
pose of confronting the VPSIC rendition. 

In FigE] the orbital-resolved DOS of YTO is showed. 
The occupied DOS have two major contributions: at 
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FIG. 5: (Color on-line) Orbital-resolved DOS for FM YTO. 
For clarity only Ti 3d, Y 5d, and O 2p are showed (oxygens 
on-top and in plane with Ti are labelled Ot and Op, respec- 
tively). Notice that Y and O DOS are magnified by more 
than one order of magnitude with respect to the dominant Ti 
3d DOS. 



VBT there is a ~0.8 eV-wide fully spin-polarized DOS 
peak of Ti 3d states (residually hybridized with a small 
O 2p portion). Despite the nominal Ti 3+ d 1 configura- 
tion, a certain amount of Ti d-0 p hybridization is clearly 
visible in Figure (notice however the different scale of Ti 
3d and O 2p DOS: here the O 2p weight is way smaller 
than, e.g. in manganites). It follows that the calcu- 
lated static charges and magnetic moment differ consid- 
erably from their nominal values (for Ti we obtain ^1.6 
and ~0.7 electrons for up and down-polarized 3d state, 
respectively, and M=0.92/xb, a bit larger than the ob- 
served 0.8/^s). Below the Ti 3d peak there is a broader, 
unpolarized DOS of O 2p character spanning the energy 
interval between -4 eV to -8 eV (not showed in Figure). 
The CBB bands are also dominated by far and large by 
Ti 3d t2g s tates, residually hybridized with O 2p and 
Y 4d orbitals. Thus we can unquestionably categorize 
the system as a true Mott-Hubbard insulator, at vari- 
ance with most manganites or cuprates that are actually 
charge-transfer insulators or in the intermediate regime 
(we will be back later on this important point). 

In the band energy plot of FM YTO (Fig® left panel) 




FIG. 6: (Color on-line) Top: Band energies at x=l/4 
(h=0.125) for (lxl) PM phase (al) and (2x2) AF phase 
(bl). K-points coordinates (units of l/a=l/b, 1/c with 
a, b, c unit-cell parameters) are X=[7r/2,0,0], X=[0,7r/2,0], 
L= [71-/2,7172,0], M= [71-/2,7172,7172]. Bottom: calculated FS for 
PM (a2) and AF (b2) phase. 



we see four occupied bands (one for each Ti) separated 
from the 3d empty conduction bands by 1.8 eV. The fun- 
damental gap only involves bands of majority manifold, 
and is direct at T. The CBB bands span a ~ 1 eV wide 
energy interval. According to our calculations, the sharp 
DOS peak at the valence top is a complex admixture 
of the five Ti 3d orbitals. To rationalize quantitatively 
the identity of this state, we have diagonalized the corre- 
sponding P^ m / density matrix in the 3d orbital subspace. 
The results are reported in Tab lllll for two coordinate 
systems: the orthorhombic Pnma \/2 x \/2x2 (x',y',z'), 
and the conventional cubic(x,y,z), which differ by a 45° 
rotation 68 of the(x,y) plane. Focusing first on the Ti sited 
at (0,0,0) in the cubic reference system of YTO, we see 
that it shows a prevailing contributions of \yz) and (\xy)) 
orbitals; however, not completely discardable e g contri- 
butions are present as well. The charge density isosurface 
plot (FigJS]left panel) shows that this state can be sub- 
stantially expressed as ~ 0.75 \yz) + 0.56|a;y). Also 
evident is the resulting orbital ordering: co-planar states 
shows an altcrnance of dominant \yz) and \xz) contri- 
butions, plus a change of sing for \xy) which makes the 
lobes of \yz) (or \xz)) leaning back and forth towards the 
(x,y) plane (e.g. |* 2 ) ~ 0.75 \xz) - 0.56|xy)). On the 
other hand, states aligned along z only differ by the al- 
ternance of \xy) sign, thus ^3) ~ 0.75 \yz) - 0.56\xy), 
1^4) ~ 0.75 \xz) + 0.56|xy). Our results are in excellent 
agreement with the finding of linear dichroism x-Ray ab- 
sorption measurement^ which gives 0.8 and 0.6 for the 
coefficients of the two most occupied ta g orbitalsJ^, and 
with LDA+DMFT resultslH which obtain 0.78 and 0.62, 
respectively. 

Now we move to analyze the results for LTO; remark- 
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FIG. 7: (Color on-line) Orbital-resolved DOS for AF G-type 
LTO (right). Orbital labels are the same as in Fig(5] La and 
Ot states are spin-compensated, due to AFG symmetry. La 
and O DOS are magnified by more than one order of magni- 
tude with respect to the dominant Ti 3d DOS. 



TABLE III: 3d orbital decomosition of the four occupied 

states (one for each Ti) at VBT of YTO and LTO. Coor- 
dinates (x',y',z') and (x,y,z) refers to orthorhombic and con- 
ventional cubic cartesian axes, respectively, as indicated in 
Figi] 

\x'y') \x'z') \y'z>) \z>*) \x' 2 -y'i) 

YTO 

Ti 1 0.11 0.48 0.58 0.33 0.56 

Ti 2 0.11 0.48 -0.58 -0.33 -0.56 

Ti 3 -0.11 0.48 0.58 -0.33 -0.56 

Ti 4 -0.11 0.48 -0.58 0.33 0.56 
LTO 

Ti 1 0.02 0.15 0.78 0.08 0.60 

Ti 2 0.02 0.15 -0.78 -0.08 -0.60 

Ti 3 -0.02 0.15 0.78 -0.08 -0.60 

Ti 4 -0.02 0.15 -0.78 0.08 0.60 

~ \xy) \xz) \yz) \z A ) jg^g) 
YTO 

Ti 1 0.56 -0.07 0.75 0.33 0.11 

Ti 2 -0.56 0.75 -0.07 -0.33 0.11 

Ti 3 -0.56 -0.07 0.75 -0.33 -0.11 

Ti 4 0.56 0.75 -0.07 0.33 -0.11 
LTO 

Ti 1 0.60 -0.45 0.66 0.08 0.02 

Ti 2 -0.60 0.66 -0.45 -0.08 0.02 

Ti 3 -0.60 -0.45 0.66 -0.08 -0.02 

Ti 4 0.60 0.66 -0.45 0.08 -0.02 




FIG. 8: (Color on-line) Charge density isosurface n±=± 0.01 
electrons/cm 3 ) of the upmost occupied state for FM YTO 
(left) and AF-G LTO (right). Red (light) and blue (dark) 
surfaces represent spin majority (+) and minority (-) contri- 
butions, respectively. On this scale only Ti d contributions are 
visible (oxygen contributes residually, see the DOS in Figs[5] 
and 17). Both YTO and LTO are orbital-ordered, i.e. the four 
Ti atoms in the cell have same integrated charge but differ- 
ent orbital distribution (numbers connect each Ti with the 
corresponding 3d orbital decomposition reported in Tab llllf) . 



able differences from YTO emerge from the calculated 
DOS (FigJT]) and band structure ©: the fundamental 
band gap is a bit smaller for LTO but still quite sizeable 
(1.6 eV); furthermore, the latter presents a band flatness 
which is stronger than in YTO: the occupied 3d states 
at VBT now span a much narrower energy range (0.4 eV 
instead of 0.8 eV), and the hybridization with the oxy- 
gens is more residual, although still well visible. Even 
the conductions bands in LTO appear flatter, and in fact 
they are separated in two groups by a gap of 0.2 eV. The 
magnetic moment is 0.89 fj,B, similar to YTO. 

The difference with YTO is also well spelled by the 
analysis of the diagonalized density matrix in Table IIIII 
Looking at the cubic reference system, at variance with 
YTO we see that the e 9 contribution is now almost van- 
ishing, and the occupied states are almost purely t2 9 . 
Moreover, the diversification of the t2 9 occupancies is 
much reduced with respect to YTO as a results of the 
smaller rotations, and indeed this state approximately 
resembles cygar-shaped [lll]-directed lobes resulting by 
the near even t2 9 combination, as confirmed by the cor- 
responding charge density isosurface plot in FigjSJ right 
panel (notice that if t2 9 coefficients were exactly the 
same, and ^2, as well as ^3 and ^4, would have been 
identical, and the resulting cygars in each plane parallel 
to each other, pointing all along [111]). 

While the orbital charge distribution in YTO and LTO 
is so different, and causes much of their macroscopic dif- 
ferences, the relative ordering is the same: even for LTO, 
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FIG. 9: (Color on-line) Pnma structure of YTO (left) and 
LTO (right). Cell parameters are fixed to experimental val- 
ues, while atomic positions were relaxed according to VPSIC. 
Labels indicate Ti-O-Ti angles and Ti-0 distances in plane 
(dp, dp) and along z (9 Z , d z ). Results are reported in Tab I IV I 



TABLE IV: Atomic positions in crystal coordinates (x/a, x/b, 
x/c), and main structural parameters (Ti-O-Ti angles in plane 
(Op) and along z (8 Z ), Ti-0 distances along z (dz) and in-plane 
(the shorter ds and the longer dz,) bonds) for Pnma YTi03 
and LaTi03 calculated by VPSIC, in comparison with the 
experimental data (in parenthesys) . Cell structures are fixed 
to the experimental values a=5.316 A, 6=5.679 A, c=7.611 A 
for YTO, a=5.640 A, 6=5.584 A, c=7.896 A for LTO^. 





x/a 


y/b 


z/c 


Y 


0.478 (0.479) 


0.073 (0.073) 


1/4 


Ti 








1/2 


Oj 


-0.139 (-0.121) 


-0.063 (-0.042) 


1/4 


On 


0.307 (0.309) 


0.184 (0.190) 


0.067 (0.058) 


La 


0.491(0.493) 


0.053(0.043) 


1/4 


Ti 








1/2 


Oj 


-0.080(-0.081) 


-0.008(-0.007) 


1/4 


On 


0.0288(0.291) 


0.204(0.206) 


0.042 (0.043) 




ds 


dz, 


d z 


YTO 


2.0(2.02) 


2.13(2.08) 


2.07(2.02) 


LTO 


2.02(2.03) 


2.06(2.05) 


2.02(2.03) 




dp 


i 





YTO 140.41° (143.62°) 133.30° (140.35°) 
LTO 153.82° (152.93°) 154.30° (153.75°) 



in the plane there is perfect alternance (i.e. chessboard- 
like order) of leading \xz) and \yz) contributions (this is 
less evident than in YTO since the t2 9 coefficients are not 
as different as in YTO), plus a sing alternance for \xy). 
Along z only the sing alternance occurs. Our calculated 
%2g coefficients are remarkably close to the values (0.56, 
0.45, 0.69) measured by NMR spectra in ReflhSl. as well 
as those calculated by a model Hamiltonian (0.6, 0.39, 
0.69) in ReflH. 

The observed magnetic ground-state is correctly pre- 
dicted for both systems: for YTO we found the FM en- 
ergy lower than AF-G and AF-C phases by 10.1 meV/f.u. 
and 8.3 meV/f.u., respectively (in agreement with pre- 
vious LDA+U results 31 with U-J=3.2 eV). For LTO, 
on the other hand, we obtain the AF-G phase lower 
than FM and AF-A phases by 15.2 meV/f.u. and 10.05 
meV/f.u., respectively. Fitting the energies on a 2- 
parameter nearest-neighbor Hciscnberg Hamiltonian: 

H = — — jVp/^i • Si+ X + S{ ■ Si+y) + J z Si ■ Si+ Z ~^20) 

i 

where i+x, i + y, and i + z indicate nearest neighbors of 
i in the x,y, and z directions, respectively, we obtain 3 p i = 
4.15 meV and Jz=1.8 meV for planar and orthoghonal 
exchange interaction parameters in YTO, respectively; 
Jpi= -5.02 meV and J z =-5.03 meV for the same quanti- 
ties in LTO. These results nicely confirm the conjectures 
derived by the analysis of the orbital ordering: while a 
remarkable anisotropy is present in YTO, LTO is sub- 
stantially isotropic. 

Table [5] shows experimental and VPSIC-calculated 
atomic coordinates and the most important structural 
parameters, i.e. Ti-O-Ti angles (8), and Ti-0 distances, 



indicated in FiglHl In plane there are two types of Ti-0 
bonds, long (dz,) and short (ds), which alternate along 
both x and y directions (see FigJS]), while along z there 
is only one d z ~ ds- These values easily rationalize the 
chessboard-like Ti d ordering: on each Ti the occupied 
state prefers to lie along the longer Ti-0 bond (thus al- 
ternatively \xz) and \xy) for dz, parallel to x, or \yz) and 
\xy) for dz, parallel to y). For YTO the difference be- 
tween dz, and ds is quite sizeable, and give rise to a very 
pronounced ordering, as seen in the analysis of the charge 
density. For LTO the dz, and ds difference is much re- 
duced, and so is the planar chessboard ordering, indeed. 
Notice that for both materials the JT-type character of 
the structural distortions is quite residual, i.e. proper- 
ties along x, y, and z are, on average, almost the same 
(especially for LTO). The GdFeO-type tiltings and rota- 
tions, on the other hand, are quite sizable and represent 
the major factor determining the observed structures and 
the consequent splitting of the tz g triplet state. Finally, 
VPSIC-calculated structure is satisfactorily close to the 
experimental determination for both LTO and YTO (al- 
though for the latter oxygen rotations are a bit overem- 
phasized along z axis). 

In summary, our VPSIC results furnish a coherent 
guideline to understand (at least part of) the differences 
between YTO and LTO, and correctly describes the dif- 
ferent magnetic ordering of the two systems: The bigger 
GdFeO-type distortions of YTO produce crucial differ- 
ences in electronic and magnetic properties, as evidenced 
by our results: a) larger Ti 3d-0 2p and t2 ff -e ff mixing; 
b) crucially different charge density distrubution around 
Ti; c) an increase by factor ^2 of the occupied 3d state 
bandwidth. The wider rotations, in particular, destibi- 
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lize the AF superexchange coupling which prevails in a 
purely d 1 t2 9 unrotated Pnma environment. 

Notice that the above considerations are completely 
reverted for doped manganites, whose chemistry is gov- 
erned by e 9 : in that case cubic symmetry and absence 
of octahedral rotations works in favor of e 9 -p hybdridiza- 
tion. In titanates, on the other hand, absence of octa- 
hedral rotations means vanishing p-d hybridization, pure 
2 g charge character, and minimal i g bandwidth. 

It remains to explain the large difference between our 
calculated and the measured energy gap. This actu- 
ally occurs by construction: our VBT and CBB band 
energies represent removal and additional energies, and 
their difference estimates the on-site Coulomb energy U, 
whereas the lowest excitation measured for these true 
Mott-Hubbard insulators is an intra-site excitation en- 
ergy which of course does not include U. The attempt 
to argument a presumed smallness of U on the basis of 
the very tiny energy gap of LaTi03 23 is a misinterpre- 
tation. In fact, according to our band structure U ^3 
eV, as expected for a system of this kind. It is also not 
very proper the strategy carried out in several works of 
estimating the excitation energy as LDA-calculated ti g 
(average) band splitting: this is justified by the fact that 
in the limit of vanishing U (i.e. delocalized electrons) 
excitation and additional/removal energies go back to be 
the same quantities; however we must keep in mind that 
here the vanishing of U is an artifact of the LDA, not 
a true feature of the titanate. A more rigorous strategy 
to evaluate the lowest excitation energy is suggested in 
refH3 where a suited Hamiltonian for the excited state is 
constructed in such a way to project out the electronic 
ground state. Here we do not pursue this route which 
overcomes the capability our actual methodological set- 
ting, and leave to future developments the investigation 
of the crystal-field splitting and orbital-liquid state in 
LTO by VPSIC. However, we emphasize that the ratio- 
nalization of the FM vs. AF G-type competition is cor- 
rectly described already at the level of our ground-state 
results. 



C. Magnetic manganites: CaMnOa 

CaMn0 3 (CMO) is a prototype AF G-type insula- 
tor. The nominal Mn 4+ d 3 configuration triggers AF 
superexchange coupling via the fully polarized major- 
ity t2 g orbitals, and AF semi-covalent exchange inter- 
actions through the empty e g states. The t2 9 spherical 
charge distribution favours a robust centrosymmetric oc- 
tahedral structure, and the near complete absence of ro- 
tations leaves the systems substantially cubic (a small 
Pnma distortion is actually observed, but it will not 
be considered here, since it is immaterial for magnetic 
and electronic properties). According to Goodenough- 
Kanamori rules^i, spin coupling is expected to be AF 
and isotropic (G-type). While this is indeed verified by 
a series of experiments and calculations, the determina- 
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FIG. 10: (Color on-line) Density of states of the most impor- 
tant orbitals (Mn d and O p) of AF-G CMO calculated within 
LDA, VPSIC, and VPSICq. Light (red) shadowed areas are 
for Mn t2g and O p orthoghonal orbitals; solid black lines for 
Mn e g and O p ligand orbitals. 



tion of the electronic and magnetic properties in detail 
is less clear and some discrepancies between the inter- 
pretation of photoemission data and the band energies 
obtained by standard local functionals make the system 
an ideal case for testing new methods. Indeed, from 
the theoretical side CMO was studied in the past using 
LDA^-iS, GGAi 2 ., LDA+U^, GGA+Ui 3 ., unrestricted 
Hartree-Fock (HF) 35 i 41 , configuration interaction (CI) 41 , 
and model (Hubbard) Hamiltonian 3 ^, while experimen- 
tally a number of optica l 24 ' 44 ' 45 i 47 and transport^ mea- 
surements has been carried out. 

The general characteristics can be illustrated with the 
help of our calculated DOS (Fig fTU)) and band energies 
(FigfTT]) for the observed AF G-type phase at experimen- 
tal lattice parameter. All theories (LDA, VPSIC, and 
VPSICo) describe the system as insulator, with a ^7 eV 
wide valence band manifold of mixed O p and major- 
ity Mn t2 S states. Very importantly, at variance with 
the nominal configuration, a consistent amout of filled e g 
states is present, in fact, in all the three calculated DOS. 
Below the p-d valence bands a narrow O 2s band mani- 
fold lies at about 18 eV from the VBT. Moving up in the 
energy above the fundamental gap we find distinct groups 
of majority e g , minority tig-, and minority e g states. 

Looking deeper at the DOS important differences ap- 
pear in the LSDA and VPSIC/VPSIC description: for 
both the p-d valence manifold shows a double peaked 
structure, in agreement with X-Ray (XPS) and ultravi- 
olet (UPS) photoemission 4 ^, but the orbital character of 
the peaks is different: in LSDA the big part of t2 9 spectral 
weight is right below the VBT, while the tail region from 
-5 eV to -7 eV is mainly O 2p. The VPSIC/VPSIC on 
the other hand recovers a spectral redistribution more in 
line with the observations, with prevalently O 2p states 
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FIG. 11: Band energies of AF-G type CMO in cubic FCC 
symmetry, calculated by LSDA, VPSIC, and VPSIC . K- 
points coordinates (units of 2ir/a,) are F=[0,0,0], X=[1,0,0], 
W=[l, 1/2,0], L=[l/2,l/2,l/2], K=[1,1,0]. Dominant and sec- 
ondary orbital character for each group of bands is also indi- 
cated. 



TABLE V: Equilibrium lattice parameter ao (in A) for AF-G 
and FM phases and exchange interactions J (in meV) of cubic 
CaMnOa calculated by LSDA, VPSIC, and VPSICo (experi- 
mental values are reported for comparison). J eq and J ex are 
values calculated for equilibrium and experimental ao, respec- 
tively. 







LSDA VPSIC VPSICo 


expt. 


ao 


(AF-G) 


3.75 


3.74 


3.83 


3.734 


ao 


(FM) 


3.77 


3.76 


3.88 




J ex 




-37.0 


-6.1 


+6.5 


-6.6S 


J ec 




-35.3 


-5.7 


+27.3 





near VBT, and the highest peak of t2 S states at the bot- 
tom of the valence bands. The LSDA failure is clearly re- 
lated to the insufficient t 2g spin splitting, which amounts 
to a mere 2.5 eV and leaves the majority t2 S much too 
high in the energy. In VPSIC the t% g splitting increases 
up to about 9 eV, which is consistent with the estimated 
value of U2i. Notice that a single energy parameter is not 
enough, however, to properly locate the t2 9 states, since 
a consistent portions of those is also present in the 4 eV- 
wide region below the VBT, where the p-d hybridization 
is strong. 

In LSDA we obtain an energy gap of 0.42 eV. Looking 
at the band picture, the VBT runs fiat between X and W, 
while the CBB e g is flat between T and X, in agreement 
with previous LDA calculations (see e.g|36l); in VPSIC 
the gap is 1.01 eV, and both VBT and CBB are flat be- 
tween r and X. Above the energy gap, LSDA describes 
the 2 eV-wide majority e g bands overlapped with the very 
narrow minority t 2g peak, whereas in VPSIC/ VPSICo 




lattice constant (A) 

FIG. 12: (Color on-line) Total energies per cell vs. lattice 
parameter calculated by LSDA, VPSIC, and VPSICo for AF- 
G (solid symbols) and FM (open symbols) ordering. Notice 
that differences in total energy due to different methods have 
no physical meaning and have been arbitrarily translated in 
the Figure for the sake of clarity. 



the latter lies about 2 eV above the centroid of the e ff . 
While we could not find in literature a clear determi- 
nation of the band gap value, interband transitions ex- 
tracted from photoemissioni^ and optical conductivity 
measurements^ seem to be very consistent with the VP- 
SIC calculation. Specifically, the distance between O 2p 
and majority e g peaks (~3 eV) and between O 2p and 
minority e g peak (~ 6.5 eV) compares excellently with 
the respective values 3.07 eV and 6.49 eV extracted by 
Lorentz oscillator fitting of conductivity spectra reported 
in Ref,45. Also quite consistent, albeit with a bit larger 
value (3.7 eV) for the O 2p-Mn e Q transition, is the va- 
lence band spectra deduced in Ref|44| by fitting a CI clus- 
ter model to XPS data. 

Now we move to examine structural and magnetic 
properties. In FigU^] total energies as a functon of lat- 
tice parameter for AF-G and FM phases are reported, 
calculated with LSDA, VPSIC, and VPSIC . Values of 
the equilibrium lattice are reported in TabfVl The trend 
is similar to that seen in Section IIII Al for wide-gap ox- 
ides: VPSIC and VPSICo respectively reduces and ex- 
pands the volume with respect to the LSDA value. Here 
however, the correction of VPSIC is very tiny, and both 
LSDA and VPSIC gives lattice constant in very good ac- 
cord (within 1%) with the experimental value. On the 
other hand, the VPSICo is much less satisfying, and the 
predicted equilibrum lattice overestimate the experiment 
by ~2.5%. For what concern the difference between AF- 
G and FM energies, LSDA is known to exaggerate the 
contribution of AF superexchange due to the excessive 
t 2g -p hybridization in the region near VBT (discussed 
in FiglTUJ) , thu predicting a strong AF-G stability (in 
agreement with previous LSDA calculations^^ through 
the whole examined range of lattice values. The VPSIC, 
on the other hand, describes a much tighter competi- 
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FIG. 13: (Color on-line) Exchange-interaction parameter J 
calculated in LSDA, VPSIC, and VPSICo; see text for the 
exact definition of J. 

tion: while at equilibrium the AF-G phase is stable, a 
moderate lattice stretching (about 1% tensile strain) is 
sufficient to reverse the ordering and stabilize the FM 
ground state. Finally, the VPSICo apparently performs 
very poorly for magnetic coupling: it completely reverses 
the LSDA description, and predicts the FM phase as ro- 
bustly stable in a wide lattice range. 

The magnetic competition can be better appreciated 
in terms of exchange-interaction parameters J, plotted in 
FigHSl (for the sake of comparison we have adopted the 
same definition of J given in ReffUl. based on the single- 
parameter Heisenberg Hamiltonian H = — Jj^y ^» ' 
where is the versor of the i-site spin direction). Val- 
ues calculated at equilibrium (J eq ) and experimental 
(J ex) structure are reported in Table |V] We can re- 
mark the excellent agreement of the VPSIC value (-6.1 
meV) with J=-6.6 meV extracted from the diagrammatic 
Rushbrook-Wood formula^ for the magnetic susceptibil- 
ity corresponding to the experimental Neel temperature 
Tat=130 K. The VPSIC also compares fairly well with 
those calculated by CI (8.1 meV) and HF (10.7 meV) 
in RefEH. On the other hand, both LSDA and VPSIC 
largely deviates from these estimates, albeit in opposite 
directions. 

It is interesting speculating on the remarkably differ- 
ent magnetic behaviour described by the three methods. 
This mainly reflects the difference in the %2g spectral 
weight seen in FigfTUl specifically the substantial DOS 
shift from the top to the bottom of the p-d band manifold 
when moving from LSDA to VPSIC and to VPSIC . As 
a quantification of this effect, in Figfrj]we report the in- 
tegrated charge of majority t^g and e g bands as described 
by the three methods: all methods describe a filled (i.e. 
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FIG. 14: (Color on-line) Integrated DOS of t2 9 (green) and 
e g (black) majority orbitals, normalized to their respective 
degeneracies 3 and 2. The DOS are also reported for clar- 
ity. Lenght of red arrows indicate the amount of t2 9 charge 
without the lowest t2 9 peak. 



integrated to 3) t2 g state, but if we evaluate the amount 
of t% g charge located in the upper part of the valence 
band region which does not include the towering lower- 
end peak and mainly hybridizes with oxygens (quanti- 
fied in Fig03] by the vertical red arrows) we find that 
this charge in LSDA is about 73%, in VPSIC just 43%, 
and finally in VPSICo a mere 33% of the whole majority 
t2 g charge. Our interpretation is the following: the lower- 
end t2 g peak is too low in energy and too little hybridized 
with oxygens to give a meaningful superexchange contri- 
bution (through t2 g - 1 poR 7r-type bonding). Thus while 
the t2 g spectral weight in transferred to the lower end of 
the valence band, the e g charge distribution remains sub- 
stantially similar across the three renditions, and its FM 
contribution takes place and even becomes dominant ac- 
cording to VPSICo- Notice that the e g charges integrates 
to a remarkable 1/2 electron per orbital at the VBT (in 
agreement with previous calculations^) thus it represent 
a prototypical case of FM coupling (via e g -pLi hybridiza- 
tion) according to Goodenough-Kanamori rules^i. 

In summary, the VPSIC seems to furnish a very con- 
sistent description of CMO for what concern electronic, 
structural, and magnetic properties. Values calculated 
in VPSIC compares well with the experiments and with 
results of other beyond-local approaches, when available, 
correcting the most obvious deficiencies of LSDA. On the 
other hand the VPSICo, while furnishing a band spec- 
trum substantially similar to that of VPSIC, fails in the 
precise account of structural and magnetic properties. 
In fact, our study revealed that subtle differences in elec- 
tronic properties can result in visible errors in some ob- 
servable quantities: specifically a slightly narrower p-d 
valence bandwidth and an excessive t% g localization to- 
wards the lower end of the valence band manifold can 
result in a 2-3% overestimation of the lattice constant 
and in a unreliable value of the exchange-interaction pa- 
rameter. 
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TABLE VI: Calculated bond-lengths of selected bonds from several molecules compared to experimental values. Both VPSIC 
and LSDA bond- lengths are shown. 8b l in column 6 (7) denotes the absolute percentage difference between the calculated 
VPSIC (LSDA) and experimental bond-length 



Molecule 


Bond 




Bond-length (A) 




5rt (%) 




VPSIC 


Experiment 


LSDA 


VPSIC 


LSDA 




B-Cl 


1.748 


1.742 


1.754 


0.349 


0.706 


CH4 


C-H 


1.081 


1.087 


1.121 


0.567 


3.090 


CofWEthanp") 


C-C 


1.478 


1.536 


1.524 


3.751 


0.776 


v .. _j. ± \ ^ y VjiuuLii cx±i ^ / 


C-C 


1.503 


1.555 


1.548 


3.315 


0.419 




C-C 


1.471 


1.501 


1.519 


1.986 


1.168 


^3 


O-O 


1.224 


1.278 


1.273 


4.200 


0.364 


NaCl 


Na-Cl 


2.317 


2.361 


2.336 


1.883 


1.039 


S1H4 


Si-H 


1.431 


1.480 


1.518 


3.344 


2.594 


SiCk 


Si-Cl 


2.020 


2.019 


2.054 


0.056 


1.747 


PH 3 


P-H 


1.407 


1.421 


1.460 


0.997 


2.736 


PF 3 


P-F 


1.576 


1.561 


1.644 


0.957 


5.349 


SH 2 


S-H 


1.341 


1.328 


1.382 


0.984 


4.030 


CuF 


Cu-F 


1.763 


1.745 


1.735 


1.009 


0.548 


ZnH 


Zn-H 


1.527 


1.595 


1.629 


4.260 


2.130 


C2H4 (Ethylene) 


C=C 


1.309 


1.339 


1.351 


2.237 


0.903 


CO 


c=o 


1.124 


1.128 


1.153 


0.384 


2.217 


CO2 


c=o 


1.142 


1.162 


1.185 


1.706 


1.959 


2 


0=0 


1.188 


1.210 


1.227 


1.793 


1.388 


N 2 


N=N 


1.086 


1.098 


1.119 


1.077 


1.924 


C 2 H 2 (Acetylene) 


C=C 


1.196 


1.203 


1.234 


0.566 


2.605 


CeHg (Benzene) 


C:C 


1.371 


1.397 


1.411 


1.871 


0.985 



IV. RESULTS: MOLECULES 

While the PSIC formulation was originally conceived 
to work for extended systems (i.e. within periodic bound- 
ary conditions) it may be no less useful for finite systems 
as well. Indeed, if for isolated atoms the full PZ-SIC is 
easy and straightforward, for large molecules and clus- 
ters its application may become methodologically cum- 
bersome and computationally extensive, and the PSIC 
can furnish a practical and reliable alternative. The im- 
plementation in local orbital basis set and pseudopoten- 
tials (ASIC), carried out in ReffHTl in the framework of 
the SIESTA code^i, is ideally suited to the aim since it 
can treat both extended and finite systems on the same 
foot. (In principle even the plane- wave implementation, 
albeit much less efficiently, can be applied to finite sys- 
tems by supercell approach.) In the last few years a 
series of works related to molecules have been carried 
out by ASIC 7 with tipically satisfying results (provided 
that the relaxation parameter a is kept fixed to unity). 
The VPSIC implemented in local orbital basis set (i.e. 
the variational generalization of the ASIC) is expected 
to yield KS spectra largely similar, although not identi- 
cal, to those obtained with ASIC. Additionally, the per- 
formance of the VPSIC energy functional (eqn. [TJ and 
the associated forces (eqn. I18[) for equilibrium molecular 
geometries needs to be investigated. 

In this section we look at the VPSIC description of the 
spectral and geometric properties of several molecules 



selected mostly but not exclusively from the standard 
G2 set 70 . Calculations were carried out using a de- 
velopment version of the SIESTA code£- within which 
the VPSIC method was implemented. Some details re- 
garding the implementation are given in the appendix. 
For all of the atomic species, standard norm-conserving 
pseudopotentials generated using the Troullier-Martins 
scheme^ are employed including core-corrections where 
necessary. Scalar relativistic pseudopotentials are used 
for the III period elements. A numerical double-zeta- 
polarized (DZP) basis selj^ is employed for all of the 
atomic species. An energy shift of 50 meV is used to 
set the cut-off radius for the pseudo atomic orbital basis 
functions. Geometry optimizations are performed using 
a conjugate gradients algoritm until all of the forces are 
smaller than 0.01 eV/A. 

A. Equilibrium bond lengths 

Table I VII shows the equilibrium bond-lengths ob- 
tained within VPSIC of selected bonds in several gas 
phase molecules. The representative set chosen includes 
molecules mainly built from I, II and III period ele- 
ments as well as non-magnetic transition metals and fur- 
thermore includes several species hosting different types 
of chemical bonds. Also presented for comparison are 
the corresponding LSDA and experimental bond-lengths. 
We find that the calculated VPSIC bond-length is gen- 
erally shorter than the corresponding LSDA estimate. 
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From column 5 in table I VII we see that the LSDA bond- 
lengths are typically a few percent longer than those 
from experiment. On the other hand (see column 4 
in table IVI|) VPSIC bond-lengths are seen to be a few 
percent shorter relative to experiment. In columns 6 
and 7 we show the absolute percentage error 5 l BL (X) = 
\l % (x) |xioo - n calculated bond-lenghts L^X) 

relative to experiment L^ xpt for each molecular species i 
and X e {VPSIC,LSDA}. The estimated mean absolute 

percentage error over the test set Asi(X) = ^ i ^ 1 ^^ — - 
comes out to be 1.84% in LSDA and 1.77% in VPSIC. We 
also note further that within the test set, the maximum 
percentage error observed within VPSIC is ~ 4% for the 
case of ZnH and that for the majority of molecules the 
error is typically under 2.0% 



B. Ionization potentials 

As in the case of the ASIC method, the primary ad- 
vatage over LSDA afforded by VPSIC is expected to lie 
in the systematic improvement of KS eigenvalue spec- 
tra. The method is thus particularly relevant for DFT 
based electron transport schemes wherein an accurate de- 
scription of the KS spectr a 73 ' 74 is often important. In 
exact KS DFT only the highest occupied orbital eigen- 
value ( e HOMO ) has a rigorous physical interpretation 
and corresponds to the negative of the first ionization 
potentia l 75 ' 76 . In general, for a N electron system, the 
following equations hold in exact KS-DFT 



HOMO 



(M) = -I N for (N — 1 < M < N) (21) 



HOMO 



(M) = -A N for (N < M < N + 1) (22) 



where —In and —An are the ionization potential (IP) 
and the electron affinity (EA) respectively. However, 
semilocal approximate functionals such as LSDA/GGA 
are known to perform poorly with regards to satisfying 
equations |2"T1 and |2"21 especially for molecular systems. In 
the following, we assess the mapping between electron re- 
moval or addition energies and the KS spectrum obtained 
from VPSIC also showing the corresponding LSDA re- 
sults for comparison. Furthermore, for both LSDA and 
VPSIC, the molecular geometries used to estimate the IP 
and EA are the corresponding equilibrium geometries in 
the neutral configuration. 

In table I VIII and figure [TS] we compare the experi- 
mental IP for several molecules with the correspond- 
ing negative e HOMO obtained using LSDA and VPSIC. 
It is clear that LSDA underestimates the removal en- 
ergies significantly in all the cases. In contrast, the 
mapping between the experimental IP and _ e HOMO from 
VPSIC is excellent. Indeed, the mean absolute devia- 
tion A IP (X) (X = LSDA, VPSIC) from experiment, 



TABLE VII: Experimental Ionization potential (IP) com- 
pared to calculated negative HOMO eigenvalues for neutral 
molecules. Columns 2 and 3 present the results from LSDA 
and VPSIC respectively. The experimental data are taken 
from reference^. 



Molecule 
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7 OQ 


1 1 fiQ 

11. Dy 


1 n ar\ 
lU.oU 




7 t<R 
1 .00 


1 *3 ^7 
lo. O t 


1 9 7^ 
iz. / 


NaCl 


5.04 


8.85 


9.80 


S1H4 


8.50 


13.75 


12.30 


SiCl 4 


7.97 


12.46 


12.06 


PH 3 


6.84 


10.96 


10.59 


PF 3 


8.32 


12.96 


12.30 


SH 2 


6.09 


10.77 


10.50 


CuF 


5.53 


11.56 


10.90 


C2H4 (Ethylene) 


6.85 


10.91 


10.68 


CO 


8.80 


13.91 


14.01 


C0 2 


9.15 


15.15 


13.79 


2 


6.74 


13.55 


12.30 


N 2 


10.13 


15.42 


15.58 


C 2 H 2 (Acetylene) 


7.16 


11.31 


11.49 


CeHe (Benzene) 


6.63 


10.51 


9.23 



Ajp(Jf) 



! (*)+ipL P ti 



N 



is estimated to be 



4.29 eV for LSDA and 0.72 eV for VPSIC (N is the total 
number of molecules in table IVIl|) . For comparison, we 
have also included in figure [T5] results obtained with a 
fully self-consistent PZ-SIC approach 77 . Somewhat sur- 
prisingly the VPSIC approximation seems to produce 
better overall agreement with experiments than the full 
PZ-SIC scheme, which is seen to overcorrect the energy 
levels. This is a rather general feature of the PZ-SIC 
scheme and it has been suggested that some additional 
re-scaling procedure is neede d 78 ' 79 . 



C. Electron affinities and the HOMO-LUMO gap 

In Hartree Fock theory, Koopmans' theorem 8 ^ implies 
that the lowest unoccupied molecular orbital (LUMO) 
energy (e LUMO ), corresponds to the EA of the N elec- 
tron system when electronic relaxation is neglected. No 
such physical interpretation exists for the Kohn-Sham 
(e LUMO ) in DFT and so the EA is not directly accessi- 
ble from the ground state spectrum of the N electron 
system. However, as equation (|2"21 indicates, the EA is 
in principle accessible from the ground state spectrum of 
the N + 1 — f (0 < / < 1) electron system and further- 
more, it must be relaxation free through non-integer oc- 
cupation. Approximate functionals such as LSDA/GGA 
however perform rather poorly even in this regard as the 
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TABLE VIII: Calculated HOMO eigenvalues for singly negatively charged molecules compared to experimental negative electron 
affinities (-EA). Columns 6,7 and 8 present the LUMO eigenvalues for the corresponding neutral species. Experimental values 
are taken from — 

Molecule eggf°(eV) ~ Exp. -EA (eV) ^ UMO (eV) ~ 





LSDA 


VPSIC 




LSDA 


VPSIC 


nc=c~ 


1.79 


-2.60 


-2.97 


-6.54 


-7.59 


CH 2 =CHT 


4.13 


-0.19 


-0.67 


-3.49 


-5.07 


HC^CO" 


2.09 


-2.39 


-2.34 


-5.98 


-8.08 


CH" 


3.80 


-0.82 


-1.24 


-4.81 


-7.15 


CH 3 " 


4.73 


-0.30 


-0.08 


-3.39 


-3.41 


CH 3 


3.52 


-1.84 


-1.57 


-5.53 


-6.15 


CH 3 S- 


2.68 


-1.37 


-1.87 


0.3 


-5.73 


HC=0~ 


5.22 


1.05 


-0.31 


-3.36 


-6.86 


CN" 


1.42 


-2.97 


-3.86 


-7.83 


-10.82 


CNO" 


1.32 


-3.58 


-3.61 


-0.62 


-4.16 


NHj 


4.54 


-0.77 


-0.77 


-4.98 


-4.33 


N0 2 " 


4.64 


-0.93 


-2.27 


-5.02 


-9.47 


OF~ 


4.95 


-2.22 


-2.27 


-2.16 


-6.14 


OH" 


4.46 


-1.52 


-1.83 


0.44 


-1.86 


PH 2 " 


2.94 


-0.98 


-1.27 


-4.55 


-4.99 


s 2 " 


2.83 


-0.01 


-1.67 


-4.5 


-6.67 


SH" 


2.54 


-1.65 


-2.31 


-0.17 


-2.72 


SiH 3 - 


2.72 


-0.60 


-1.41 


-3.85 


-5.38 




FIG. 15: Experimental negative ionization potential IP com- 
pared to the calculated HOMO eigenvalues for molecules. The 
experimental data are from reference^, while the star symbol 
represents full PZ-SIC calculations from reference—. 



N + 1 electron state is often unbound with a positive 
eigenvalue. Therefore in practice, electron affinities are 
usually extracted either from more accurate total energy 
differences 82 , or by extrapolating them from LSDA cal- 
culations for the N electron system 8 ^. The failure of ap- 
proximate functionals in reproducing the spectra of an- 
ions has been traced in most part to the SI error and 
so SIC schemes are expected to perform better in this 
regard. 



FIG. 16: Experimental negative electron affinities (-EA) com- 
pared to calculated HOMO eigenvalues of negative radicals. 



In table I VI III we compare HOMO energies (denoted as 
e N+i ) °f several singly negatively charged molecules 
with the experimental electron affinities of the corre- 
sponding neutral species. We also report the LUMO en- 
ergies for the molecules, most of which are radicals, in 
their neutral ground state (denoted as e^™ )- Relaxed 
geometries of the neutral molecule are used for both the 
neutral and charged cases. We find that various — c^+i 40 
obtained from VPSIC once again map quite well onto the 
corresponding experimental electron affinities in contrast 
to LSDA which yields unbounded states with positive 
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TABLE IX: Calculated HOMO-LUMO gaps (Eg) of neutral 
molecules from LSDA and VPSIC. <5Eg in colum 4 represents 
the difference between corresponding VPSIC and LSDA gaps. 



Molecule 


Eg 






LSDA VPSIC 


SE G 



1.92 
0.84 
1.66 
1.41 
2.61 
3.16 
2.79 
2.7 
5.43 
2.16 
2.23 
2.3 
3.85 
3.39 
1.56 
2.85 
1.84 
2.38 
1.9 
3.37 
2.5 



eigenvalues are generally significantly lower than the cor- 
responding LSDA ones (see table EEE}, the HOMO- 
LUMO gaps differ by a smaller extent. This is because 
in contrast to other methods such as LSDA+U^ wherein 
the occupied levels are pushed lower and the unoccupied 
ones are pushed higher relative to the LSDA spectrum, 
within VPSIC the whole spectrum is corrected lower with 
the occupied and empty levels being shifted by different 
amounts depending upon their orbital character. For in- 
stance, the mean absolute difference between the LSDA 
and VPSIC HOMO-LUMO gaps for the test set in ta- 
ble [IX] comes out to be ~2.52 eV while the correction to 
the HOMO levels alone with respect to LSDA is around 
4 eV. In general, VPSIC is expected to open the HOMO- 
LUMO gap substantially in systems where the occupied 
and un-occupicd KS eigenstates have markedly different 
atomic orbital projections. Finally, it is worth noting 
that in contrast to explicity orbital dependent methods 
such as PZ-SIC, VPSIC does not exhibit the derivative 
discontinuity at integer occupations. Thus the eigen- 
value of the highest occupied orbital relaxes continuously 
across fractional occupations although the range of eigen- 
value relaxation is generally smaller than in LSDA. 



V. CONCLUSIONS 



BC1 3 


4.84 


6.76 


CeHe (Benzene) 


5.33 


6.17 


C 2 H 2 (Acetylene) 


6.61 


8.27 


C 2 H 4 (Ethylene) 


5.58 


6.99 


C 2 H 6 (Ethane) 


9.03 


11.64 


CH 4 


10.63 


13.79 


CO 


6.62 


9.41 


C0 2 


8.33 


11.03 


CuF 


1.5 


6.93 


C4Hg(Cyclobutane) 


8.13 


10.29 


C3H6 (Cyclopropane) 


8.15 


10.38 


N 2 


7.97 


10.27 


NaCl 


2.9 


6.75 


2 


2.16 


5.55 


3 


1.67 


3.23 


PF 3 


6.26 


9.11 


PH 3 


6.63 


8.47 


SH 2 


5.62 


8.0 


SiCl 4 


5.89 


7.79 


SiH 4 


8.64 


12.01 


Ti0 2 


1.51 


4.01 



e^OMO £ or a ji gygtems considered. Over the set of 
molecules in table IVIII[ the mean absolute error with 
respect to experiment for the electron affinities 



&ea(X) 



E 



N 



HOMO,i 



i=l \ C N+1 



(A)+EA^ xpt | 



N 



(23) 



(X = LSDA, VPSIC), stands at 4.67 eV and 0.54 eV for 
LSDA and VPSIC respectively. In figure [TC] we present 
our data together with c^+i 10 as calculated using the 
PZ-SICffi. For EAs as well, we see that PZ-SIC seems to 
systematically overcorrect the LSDA shortfall. 

We now discuss briefly the HOMO-LUMO gap in VP- 
SIC. As is apparent from columns 5 and 6 in table IVHTl 
the LUMO eigenvalues of the neutral species differ sub- 
stantially from the corresponding negative electron affini- 
ties both within LSDA and VPSIC. In general DFT 
LUMO states are expected to be lower than -EA by an 
amount equal to the derivative discontinuity A xc defined 
as 



A xc = lim e H° MO 



-N+f 



LUMO 

-N 



(24) 



i.e. A xc is the discontinuity in the eigenvalue of the 
LUMO state at AT. Thus the HOMO-LUMO gap is usu- 
ally underestimated with respect to the true quasipar- 
ticle gap Eg = 1^ — . In table |IX] we compare the 
HOMO-LUMO gaps from LSDA and VPSIC calculated 
in the neutral configuration for the molecular test set 
of table EED We see that although the VPSIC HOMO 



In summary, we have introduced the first-principles 
VPSIC approach, a variational generalization of the 
method formerly known as P SIC/ ASIC, based on the 
idea of removing the spurious self-interaction from the 
local-density functional energy. In VPSIC the self- 
interaction is removed in effective albeit approximate (i.e. 
orbitally-averaged) manner, which gives several advan- 
tages over the full SI removal (applied e.g. in the PZ- 
SIC approach or related methods for extended systems) 
such as the conservation of translational invariance (i.e. 
Bloch theorem) and the invariance of the energy under 
unitary rotations of the occupied KS eigenfunction man- 
ifold. The VPSIC approach emerges as applicable to a 
vast series of systems (insulators and metals, magnetic 
and non-magnetic, extended or finite) with an overall 
satisfying accuracy, and furthermore not significally more 
demanding than LDA/GGA from a computational view- 
point. 

We have implemented the method in two method- 
ological frameworks: plane-waves basis set and ultra- 
soft pseudopotentials, and local-orbital basis set plus 
norm-conserving pseudopotentials. The first was applied 
here to extended systems: non-magnetic oxides, mag- 
netic titanates, and magnetic manganites, the latter to a 
vast range of molecules, testing structural and electronic 
properties. Overall, the performance of the VPSIC can 
be synthesyzed as it follows: the predicted equilibrium 
structures are substantially in the same range of accuracy 
than the LDA/LSDA, tipically reducing, on average, the 
bond lengths predicted by the latter. This corresponds 
to an average underestimate of the experimental lattice 
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constant by about 1% for bulk systems, 2% for molecule 
bond lenght. On the other hand the electronic proper- 
ties are highly improved with respect to LDA/GGA, with 
band gaps of bulks, and ionization potentials and electron 
affinities for molecules tipically within 10% from the re- 
spective experimental determinations. These preliminary 
results encourage us to pursue further explorations of the 
VPSIC approach for finite and extended systems alike. 

VI. APPENDIX I - GENERALIZATION OF 
VPSIC FORMALISM TO USPP 
IMPLEMENTATION 

For the study of large-size magnetic and strong- 
correlated systems, the ultrasoft pseudopotential (USPP) 
method 62 associated to plane- wave basis set, represent 
a formidable tool which allows the use of cut-off ener- 
gies as low as 30-40 Ryd even for 'hard' transition metal 
ions such as Mn or Cu. The trade-off for this invaluable 
computational efficiency is a sensible complication of the 
VPSIC formulas presented in Section [TTJ In the following 
we furnish the VPSIC formulation adapted to the USPP 
formalism, which is our implementation of choice. For 
brevity however here we do not revisit the basic USPP 
formulation (for which we remand the readers to the orig- 
inal articles I62I). but only focus on the additional part 
specific for VPSIC. 

In the USPP approach the atomic valence charge 
is partitioned in outer ultrasoft (US) and intra-core, 
augmented (AU) contributions, of whom only the first 
changes self-consistently with the surrounding chemical 
environment. Thus the charge associated to the Bloch 
states V'^k appearing in Section [TT] only represents the 
very smooth US part. The Bloch states obey the follow- 
ing generalized orthonormality conditions: 

«kl$K'k> = *n,n> (25) 

where the overal matrix S is given by: 

S = I + \Pa,»)qab,v(Pb,v\ (26) 

Here /3 a ,v( v ) and q_ a bv are atomic projector functions 
and augmented charge, respectively, characteristics of the 
USPP formalism, and (a,&) label atomic quantum num- 
bers (l a ,m a ) (q a bp, only for l a =h)- Consistently, the 
total charge density is generalized as: 

«W = EWk|SWICk) (27) 

rikxT 

S(v) = |r><r| +Y, \P a ,v)Q abv {v){P b , v \ (28) 



where Qabvi?) are augmented atomic charge densities. 
Within this generalized framework, the VPSIC energy 
functional described in Eq[T] only includes the ultrasoft 
(US-SIC) part. In order to recover the full VPSIC en- 
ergy functional, a further augmented contribution must 
be added: 



j^VPSIC-AU _ 1 \^ ocr p a cAU /on) 
a - ~2 2^ tS a,b, V rba,v t ba,v> 

abav 

where B° b v is the matrix of Bloch state projections 
onto the beta-function basis: 

B5b, u = £ fSk (Ckl/^XAJCk) (30) 

nk 

and £ AU are SI energy associated to the augmented 
atomic charges Q ab , v (r): 

£a b> u = fdrQ a bA*)VHXc[na»(r),0}- (31) 

(In radial simmetry VHXc[n av , 0]=VHxc[n bv , 0] for 
rn a ^ m,b). Notice that we use different indices for the 
USPP projector (a, b) and for the SI projector (i, j), 
since the two basis set are conceptually and practically 
different: the latter is built on a minimal set of atomic 
orbitals, and to be physically sound it should be associ- 
ated to bound atomic states; on the other hand in or- 
der to improve the USPP transferability it is customary 
to include in the (a, b) matrix more than one state per 
angular moment (tipically the bound atomic state plusd 
one unbound state relative to some diagnostic energy ref- 
erence). This difference introduces some ambiguities in 
Eg I3T1 relative to the definition of VHxc[n a v, 0] and P a bv 
The ambiguity can be solved by associating the same 
atomic Vhxc (relative to the bound state) to all the beta 
projectors with same angular momentum l a . Another 
possibility is rewriting the VPSIC-AU energy as it fol- 
lows: 



eVpstc-au m] = * £ p[jv pJiv € au (32) 

ijua 

where: 

<£ V = E I ) E ab,u ( Pbv I fa ) (33) 

a b 

is just the augmented-only SI energy relative to the 
atomic state i at full occupancy, and can be directly cal- 
culated in the atom. The use of the simplified Eg I3"2l 
bypass altogether the explicit presence of the augmented 
charges, thus greatly simplifying the VPSIC AU energy 
functional calculation. Since our many test cases reveal 
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that EqsJ2niand[32]give indeed very similar results, we de- 
cide to adopt the latter as standard choice. Eg [32] brings 
a corresponding contribution to the VPSIC KS Equa- 
tions: 



dE VPSW-AU _ Qp^ 

l i y Bib* 3 V 

ka ijj; rnka 



(34) 



Finally, the orbital occupations defined in Eqj2] must 
be also generalized as: 



nk/ ) 



(35) 



nk 



where the ultrasoft atomic orbitals <pi U have been re- 
placed by 



\4>iv) = E S i'l>X\^i'v) 
i' v' 

with S given in Eq l26l Thus, we have: 

\S(j>iv) = \4>iv) + E \Pan)Qab,)iWbn\(l>iu) 

and 



(36) 



(37) 



VII. APPENDIX II - FORCES FORMULATION 
WITHIN PANE- WAVES AND USPP 

In case of USPP formalism the forces espression given 
in EqfTH] (or Eq[19] if we the simplified approach is con- 
sidered) must be generalized in order to include the con- 
tribution generated by the additional AU energy in [32] 



dE AU [{iP}} 



— E /nk 



ij,nko 



P- e 



39) 



In plane waves basis set the implementation of EasfTSl 
[T9l or [35] is rather straightforward except for one ingredi- 
ent which requires attention: the atomic orbital deriva- 
tive. 

The simplest case is that of atomic orbitals which re- 
main centered on the atomic positions (i.e. orbitals which 
simply translate along with their reference atom displace- 
ment). In this case the force on a given atom R„ only 
depends on the change of the orbitals sited on v : and the 
orbital derivative is easily calculated as: 
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(38) 



is the overlap matrix in the atomic orbital basis set. By 
construction this is Hermitian and also positive-definite, 
so its square root matrix can always be defined in the 
following, unique way: since (let's drop the atomic index) 
[S^ 1 / 2 , S}=0, S can be diagonalized in the subspace, 
and from its eigenvalues Xk we have S 



-1/2 
kk' 



>kk' 



A 



-1/2 



The latter is finally rotated back to the original (i,j) basis 

— 1/2 

set to obtain &. ' . 

The replacement of simple atomic-site centered cj>i U or- 
bitals with the 4>i orbitals given by Eq l36l (known as 
Lowdin orthonormalization 66 ) is required by the neces- 
sity to enforce, at any {R}, the orthonormality condi- 
tions {4>iu\(t>3ii') = ^ij5vi>' , and in turn, the correct normal- 
ization of the orbital occupation matrix defined in Eq |35l 
(upon diagonalization) < P? iu < 1 and Pf iu =N, 

with N the total number of electrons in the cell. These 
constraints are essential to the interpretation of P? iv as 
physically sound orbital occupancies. 

On the other hand, the Lowdin renormalization implies 
a remarkable complicacy in the forces formulation, since 
it is clear from Eqsl36l and 1381 that <f>i V is not simply cen- 
tered on a single atomic site, but includes contributions 
from the overlap with all other atomic orbitals 4>j^ as well 
as beta functions (3 a/1 of the cell. Since the forces formu- 
lation in the case of Lowdin-normalizcd orbitals may be 
useful even in other methodological contexts (e.g. in the 
LDA+U method, whose Hamiltonian is also written in 
terms of orbital occupancies) we dedicate the next Sec- 
tion to describe it in full detail. 



d 



0R„ 



(k + G | <p iv ) 



d 



- l (k+G).R„ (k + G |^ o)(40) 



= -i(k + G) (k + G 



(41) 



where clearly <f>i V — (pi V (r— R„), and 4>io = 4>i{r). How- 
ever, as discussed in the previous Session, orbitals <pi U 
must be replaced by cf>i U , and the forces equation gener- 
alized accordingly: 



dE 



VPSIC 



? LSD 



0R„ 



ij,nk(T 



+ C.C. 



ijti,nk(T \ J 



+ > ; /nk \ PijA" (Ckl^K^ICk) - 



E 

ijf^nka 



dK v ' 



where the first two terms account for the US contri- 
bution, and the third for the AU part. The presence 
of Lowdin-normalized orbitals brings one more sum over 
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atomic positions in the second and third term, since now 
the displacement of one single atom in R„ changes in 
principles all the orbitals, not just those sited on R„. 
The Lowdin orbital derivatives gives: 



0R V ' 



S fv',3^ nk 



dR u ' dR v 



(V>„k|^ ?v >(43) 



(where sum over repeated indices is understood). Let's 
start considering the first term: 



,9S4> 



dR v 



= (^nk\S\- 



dR„ 



(V'n 



OS 



(A 



dR u 



av'/Qabv' \Pbv' 



dR v 



Ma 



abu' 



dR v 



)qabv'{Pbv'\<t>j>^>) + 



(lpnV.\f3av>)q a bv' 



,df3, 



bu' 



8R U 



(44) 



Here the first term in square brackets select the contri- 
bution to the derivative due to the atomic orbital ^j'ij,' 
displacement, while the second term select the contribu- 
tion due to the beta functions displacement. Despite the 
apparent intricacy Eq 1441 is rather straightforward to cal- 
culate in plane waves, since all these derivatives are easily 
obtained through EqflXl 

More involved is the calculation of the second term in 
Eg 1431 which includes an uncommon square root matrix 
derivative. We can proceed as following (dropping the 
atomic indices for brevity): from S^ 1 / 2 S~ 1 ' 2 =S~ 1 we 
obtain: 



+ 



dS 
dR„ 



(47) 



whose esplicit espression is clearly similar to what al- 
ready reported in Eql44| and we can give it as under- 
stood. So, the matrix at the left side of EqJJS] is deter- 
mined. Now looking at the right side, we see that if 
commute with its derivative, the latter is easily extracted 
as: 



as- 1 / 2 1 as- 1 



dR 



2 <9R 



s- 1 ' 2 



(48) 



However, they do not generally commute (except when 
the Hermitian matrix 5 -1 / 2 is also real) and Eql48l 
does not hold. Thus, we need to solve Eg 1431 which is 
nothing but a Lyapunov matrix Equation of the form 
B=XA+AX, where the known terms are A—S^ 1 ^ 2 and 
B=dS- 1 /dR, and X=dS- 1 ' 2 /dR. The general Lya- 
punov Equation can be solved exactly, but for the specific 
values of A and B we can apply the simple strategy pro- 
posed in Ref. [65l that we repeat here for the commodity 
of the reader: we can rewrite A=CVC + , where V and 
C can be easily determined as the diagonalized matrix 
and the basis change matrix; then we can multiply both 
members of Eq. |45]by C + on the left, and C on the right: 



C+BC = C+AXC + C+XAC = VC+XC + C+XCf49) 

introducing R=C+ BC and Y=C+XC, EqUUcan be 
rewritten as R = VY +YV, which is now trivially solved, 
given the diagonal character of V: 



Y- ■ — 



Ri 



Vi 



(50) 



33 



Once Y is determined, the unkonwn X can be finally 
obtained as X = CYC + . 



8S~ 



OR 



° .£-i/a + 5-1/2 (45) 



OR 



OR 



The left side of Eq|45] can be transformed using the 
following espression: 



os- 1 _ ! os j 



(46) 



where S 1-1 can be easily obtained from S (just like 
5 -1 / 2 , as explained in the previous Section). Then we 
need to calculate the overlap matrix derivative (reintro- 
ducing atomic indices for clarity): 



dS.. 



d~R v 



dR v 



S | <pi'fi' 



VIII. APPENDIX III - VPSIC FORMALISM 
WITHIN ATOMIC ORBITALS BASIS SET 

As the VPSIC correction is based on a projection of 
the occupied KS orbital manifold onto a localized sub- 
space of atomic orbitals, the formalism naturally lends it- 
self to implementation within a localized orbital basis set 
framework. In this appendix we provide some details of 
the current VPSIC implementation within the SIESTA 71 
code. The first step in setting up the VPSIC algroithm is 
to construct a minimal set of atomic orbitals {4>i, v } & n d 
the associated projectors {7i,„} which will be used to cal- 
culate the occupation numbers (eqn. [5]) and effective SI 
energies (eqn. [3]). Within SIESTA, the functions <f>i >v are 
numerical pseudo-atomic orbitals with a finite range con- 
structed as solutions of the atomic Schrodinger equation 
with an additional confining potential at the cutoff radius 
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r<r^. The finite extent of the functions 4>i,v a ^ so ensures 
that the corresponding SIC projectors 7$ „ (eqn. [4j also 
vanish beyond the cutoff r c . In the current implemen- 
tation, the SIC potential Vhxc{Pv,u 1] in equation @] 
is obtained from a full PZ-SIC-LSDA 64 calculation for a 
free atom and imported into SIESTA via a pseudopoten- 
tial. An appropriate choice for the cutoff radius r c is then 
dictated by the requirement that the expectation value 



The expression for the VPSIC contribution to the atomic 
forces also involves only two center integrals and their 
derivatives. Setting S a ,i V = (a\4>^ v ), S ajiv = (a\(f>i tU ) 
and Ga tiL , = (a\j i>u ), 



? sic 



dE hl0 [{jj}} 



Se 



sic 



dr<j) htmi (r) Vffxck,ii(r); 1] m ; (r) (51) 



reproduces the PZ-SIC-LSDA correction of the corre- 
sponding orbital in the free atom to within a small tol- 
erance. Simultaneously, the cutoff should be reasonably 
short so as not to change the connectivity of the matrix 
elements of the PAO Hamiltonian. Therefore, in practice, 
we set the cutoff radius for the projection orbitals <f>i V on 
a given atom to be either equal to the largest among the 
cutoff radii of the PAO basis set for that particular atom 
(typically the first Q of the lowest angular momentum), 
or, if shorter, the radius at which Sefl < O.lmRy. For 
typical cutoff radii (6 to 9 Bohr) , we find that the atomic 
SIC-LSDA eigenvalues are reproduced to within 1 to 5 
rnRy for the most extended shells and to within 0.1 mRy 
for more confined shells. Thus 5ef{, c are rather well con- 
verged already for cutoff radii defined by PAO energies 
shifts— of around 20mRy. 

Using the orbitals <^„ and the projectors {"fi^}, the 
occupation numbers pfj U and the effective SI energies 
e ijav f° r the extended system can be calculated. Dif- 
ferent choices are possible for the projection operators 
that yield the occupation numbers pfj V - In our imple- 
mentation we use the so called dual projection operator 
given by 



P% v = \{\<kv)(<l>jA 



\(f>i,u)(<t>j,v\} 



(52) 



where \4>i. v ) is the dual orbital of \4>i,v) and is given by 

with S" -1 being the inverse of the overlap matrix over the 
non-orthogonal set {(/>i, v } 

Siujij. = ((/>i,v\<i>j,ii) (54) 
The dual orbitals satisfy the orthogonality relation 



(<t>i,v I <f>j 



(55) 



With this choice for the occupation number operator, 
the matrix elements of the VPSIC potential between two 
basis functions a, (3 come out as 



ySIC _ 
v a/3a — 



9 H 4^{^[H</^)(<^ 



+P°iMli,«)Ciju{~fjAP) (56) 



e si 



dp 



de SI 

2 ^^ V ° dR u PlJV dR u 



(57) 



with 



dp] 



• j i' 



dR„ 



2 OT) [Sa,ivS@ t j v + S a ,ivSpjv\ 

a/3 M 



+ P P a[ ~dR^ b ^ + ba < i "~dR u ~ 



U 1 



or,, 



and 



de SI 



dp 



/3a 



a/3 p 



dR„ 



(58) 



dR u 



where the sum is over the SIESTA basis functions a,fi 
and p% a is the density matrix given by 



^ Q =E/«k(/3|Ck)(V^ k |a) 



(60) 
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